Ion Concentration Dynamics as a Mechanism for Neuronal Bursting 
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We describe a simple conductance-based model neuron that includes intra- and extra-cellular 
ion concentration dynamics and show that this model exhibits periodic bursting. The bursting 
arises as the fast spiking behavior of the neuron is modulated by the slow oscillatory behavior 
in the ion concentration variables, and vice versa. By separating these time scales and studying 
the bifurcation structure of the neuron, we catalog several qualitatively different bursting profiles 
that are strikingly similar to those seen in experimental preparations. Our work suggests that ion 
concentration dynamics may play an important role in modulating neuronal excitability in real 
biological systems. 
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I. INTRODUCTION 

The Ho dgkin- Huxley equations [ij are of fundamental 
importance in theoretical neuroscience. These equations 
assume that the intra- and extra-cellular ion concentra- 
tions of sodium and potassium are constant. While this 
may be a reasonable assumption for the squid giant axon 
preparation (for which the equations were originally de- 
veloped), its validity in other cases is not clear. In the 
mammalian brain, for example, the neurons are much 
smaller and they are more tightly packed, resulting in 
significantly smaller intra- and extra-cellular volumes. 
Thus, typical ionic currents can have a much larger effect 
on the ion concentrations in this case. 

The effects of extracellular potassium {[K]o) accumu- 
lation on neuronal excitability have long been recognized 
and deficiencies in [K]o regulation have been im- 
plicated in various types of epilepsy (for a review, see 0) 
and spreading depression 0,H|- More recently, computa- 
tional studies have begun to clarify the role of impaired 
[K]o regulation [g-[l5j as well as other varying ion con- 
centrations [TgI, [l7|. 

In this work, we consider from a dynamical systems 
perspective the role of ion concentration dynamics in the 
generation of periodic bursting behavior. To emphasize 
the generality of our approach, we base our model on 
the standard Ho dgkin- Huxley equations. We augment 
these with additional equations that describe the dynam- 
ics of both intra- and extra-cellular sodium and potas- 
sium. The inclusion of sodium is relatively novel and 
plays a crucial role in the dynamics described here. We 
also include terms describing pumps, extracellular dif- 
fusion, and a simple glial buffering system. A different 
analysis of this system was presented in [l8| . 



* Email address: |ebarretoOgmu.edu 



t Email address: |j cressma@gmu.edu 



II. MODEL 

We begin by explicitly adopting the standard conven- 
tion that an outward membrane current is defined as be- 
ing positive fl9l|. Thus, the membrane potential V is 
given by 
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membrane 
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where Imembrane represents the sum of the various mem- 
brane currents. We aim in this work to consider a very 
simple and general model neuron. Hence we include only 
the standard Hodgkin-Huxley sodium current (with in- 
stantaneous activation), the delayed-rectifier potassium 
current, and leak current. We write the latter in terms of 
separate sodium, potassium, and chloride contributions 
[TB]. Thus, 
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where the gi [i — Na, K, CI) are maximum conductances. 
Time is measured in milliseconds, voltage in millivolts, 
and C, /, and g are measured in units per unit of mem- 
brane area, i.e., fiF/crn^, ^A/crn? , and mS/cm?, respec- 
tively. The reversal potentials Ei are given in terms of 
the instantaneous intra- and extracellular ion concentra- 
tions by Nernst equations: 



= 26.64 In 
= 26.64 In (1^) . 
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We fix Eci = -81.9386 mV. 

The extracellular potassium and intracellular sodium 
concentration dynamics are given by 
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where the concentrations are measured in milhmolar 
(niM). 7 — 4.45 x 10~^ is a unit conversion factor that 
converts the membrane currents into mM/sec (see Ap- 
pendix), and (3 — 7 is the ratio of the intracellular to 
extracellular volume The terms with tildes are the 
molar currents (mM/sec) given below, and r = lO'^ bal- 
ances the time units. 

The pump, glia, and diffusion molar currents are given 

by 
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We set the default parameter values to p = 1.25 mM/sec, 
G = 66.666 mM/sec, and e = 1.333 Hz. kbath repre- 
sents the potassium concentration in the reservoir, i.e., 
the bathing solution for a slice preparation, or the vas- 
culature in vivo. We set kbath — 4 mM for normal physi- 
ological conditions. 

The intracellular potassium and extracellular sodium 
concentrations are obtained from the following simplify- 
ing assumptions that allow us to reduce the dimension- 
ality of the system [l^ : 



[K], = 140mM4 
[Na]o = 144mM 



(18mM - [Na]^) 
~ (3{[Nal^ - 18mM) 



(6) 



The first assumption is that the sodium membrane cur- 
rent is the dominant means by which sodium is trans- 
ported across the membrane, and that during the course 
of an action potential, the transport of sodium and potas- 
sium are simply related. The second assumption is that 
the total amount of sodium is conserved. 

The remaining parameters and the equations for the 
gating variables are given in the Appendix. 



III. RESULTS 




-100 



10 20 30 40 50 

[K]„, mM 



03 



40 
30 
20 
10 




B 



SNIC 




Hopf 



Spiking 



10 20 30 

[K]„, mM 



FIG. 1: Bifurcation diagrams describing the neuron's asymp- 
totic behavior with fixed ion concentrations. In (A), [Na]i = 
lOmM. The detailed structure in the upper part of (B) is 
shown in Figure (5] 



Fixed ion concentrations 



We begin with a discussion of the dynamical structure 
of our model subject to constant values of the ion concen- 
trations. That is, we set [K]o and [Na]i to fixed prede- 
termined values, and obtain the remaining concentrations 
using Equations ([6|). We then examine the behavior of 
the neuron as given by Equations ((H)-©. 

Figure [1] shows bifurcation diagrams obtained under 
these conditions. The features in these diagrams clarify 
the neuron's bursting behavior as the ion concentrations 
undergo slow oscillations, as explained below. 

Figure [T]4. is constructed by holding [Na]i fixed at 
lOmM and plotting the asymptotic values of the mem- 
brane potential V versus several fixed values of [K]^. 
Below approximately 5.7mM, the neuron is attracted to 



a stable equilibrium (shown as a solid line) that corre- 
sponds to the resting state. At approximately [K]o — 
5.7mM, this equilibrium coalesces with a coexisting un- 
stable equilibrium (dotted line) in a saddle-node bifur- 
cation that occurs on an invariant closed curve (of infi- 
nite period) that appears simultaneously. This scenario is 
known as a SNIC bifurcation 3J]. Beyond this, for values 
of [K]o between 5.7mM and 35.2mM, a stable limit cycle 
appears, reflecting regular spiking in the neuron. This is 
depicted in the diagram by filled circles that mark the 
maximum and minimum values of the membrane voltage 
during a cycle. For increasing values of [K]o approach- 
ing [K]o = 35.2mM from below, the amplitude of this 
periodic orbit decreases and the orbit eventually merges 
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FIG. 2: (Color) Asymptotic behavior of the ion concentra- 
tions as kbath (in mM) is varied. The curve on the left de- 
notes stable equilibria and corresponds to ktath ranging from 
4.0 to 7.615. Projections of limit cycles (loops) are shown for 
kbath = 7.62, 8.0, 8.5, and 8.95 mM; the thickness of the loops 
on the right reflect small fast ion concentration fluctuations 
due to spiking behavior. The curve in the upper right, for 
kbath ranging from 8.8 to 15.0, denotes small-amplitude loops 
that correspond to tonic spiking (see inset). 



with the coexisting unstable equilibrium in a supercritical 
Hopf bifurcation. For [K]o > 35.2mM, a stable equilib- 
rium is found - this is the state of depolarization block 

Figure [1)3 is a two-dimensional bifurcation diagram 
which shows the location of the SNIC and Hopf bifur- 
cations as the value of [Na]i is varied. These curves de- 
lineate the boundaries of different attracting behaviors 
of the neuron. To the left of the SNIC curve, the neu- 
ron is attracted to the resting equilibrium. Between the 
SNIC and the Hopf curves, the neuron exhibits regular 
spiking, and to the right of the Hopf curve, the neuron is 
attracted to the depolarization block equilibrium. (The 
detailed structure at the top of this diagram is discussed 
below.) 



B. Dynamic ion concentrations 

We now describe the behavior of the full system, in 
which the ion concentrations are allowed to evolve dy- 
namically. In Figure [21 we plot the asymptotic behavior 
of the ion concentrations for several values of kbath ■ Also 
included in the figure is a portion of the SNIC curve from 
Figure [Tj3. 

At the default parameter values described above (with 



kbath = 4.0 mM), the entire system approaches a stable 
equilibrium resting state for which the membrane voltage 
and the ion concentrations assume fixed values. As kbath 
is increased, these equilibrium values change and sweep 
out the solid curve shown on the left of Figure [2j At 
approximately kbath — 7.615 mM, this curve collides with 
the SNIC boundary. Just beyond this value, the system 
jumps to a limit cycle. As kbath continues to increase, the 
projection of this limit cycle onto the ion concentration 
variables drifts upward and to the right, as shown in the 
figure. Henceforth, for brevity, we refer to such limit 
cycle projections as "loops". 

Note that these large-amplitude loops straddle the 
SNIC curve. In addition, they have periods on the order 
of several tens of seconds. Consequently, as the system al- 
ternately transitions between the resting and the spiking 
regions, there is ample time to exhibit those asymptotic 
behaviors. That is, the neuron bursts. For example, for 
kbath = 8.0 mM, the ion concentrations are attracted to 
the second (red) loop in Figure [H and the correspond- 
ing behavior of the membrane voltage is shown below in 
Figure [IK. 

At approximately kbath — 9.0 mM, the large-amplitude 
loop disappears. For larger kbath, the ion concentrations 
exhibit very small-amplitude loops, as shown in the inset 
of Figure [2] Since this loop lies entirely within the spik- 
ing region, the neuron exhibits tonic spiking. Indeed, the 
loop itself represents the small changes in the ion concen- 
trations due to individual action potentials. 

We note in passing that we have observed multista- 
bility [23l - [25j . For values of kbath approximately in the 
interval (8.8, 9.0) mM, a bursting solution coexists with 
a tonically-spiking state for the same parameter values. 
This is consistent with the analysis based on a reduced 
model in 



C. Catalog of bursting types 

The results presented above demonstrate that for pa- 
rameter values in appropriately chosen ranges, the sys- 
tem evolves on a limit cycle whose projection onto the ion 
concentration variables forms a loop that straddles the 
SNIC curve, and the neuron bursts. Bursting behaviors 
of various qualitatively different kinds can be exhibited 
by the system if similar ion concentration loops straddle 
the bifurcation curves of Figure[T}3 in different ways. Ac- 
cordingly, we can catalog all the possible arrangements, 
and examine the nature of the resulting bursting pat- 
terns. 

Figure [3| shows four different ion concentration loops 
labeled A, B, C, and D. Loop A is the kbath = 8.0mM 
loop discussed above, and the corresponding membrane 
voltage trace is shown in Figure Loop B {G = 20, 
e — 0.133, and kbath ~ 22) straddles both the SNIC and 
the Hopf bifurcation curves, and hence displays a qualita- 
tively different bursting pattern. As the loop is traversed, 
the ion concentration trajectory moves from the rest- 
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FIG. 3: (Color) Loops A-D represent the time evolution of the 
ion concentrations as the system exhibits limit cycle behav- 
ior. The loops are traversed in a counter-clockwise manner. 
The dotted and solid lines are the SNIC and Hopf bifurcation 
curves, respectively. 



ing region into the spiking region by crossing the SNIC 
curve, and then continues across the Hopf curve into the 
depolarization block region. It then bends around and 
crosses these regions in reverse order, and the cycle re- 
peats. Consequently, the membrane potential shows a 
bursting pattern that moves from quiescence to spiking 
to depolarization block and back again, as is shown in 
Figure B^. Loop C (G = 6, e = 0.7, and hath = 22) 
straddles only the Hopf curve, and hence the membrane 
potential displays round-shaped bursts, reflecting the su- 
percritical nature of the Hopf bifurcation, as shown in 
Figure dUC). Finally, Loop D (p = 0.9, G = 10, e = 0.5, 
kbath = 20, and 7 = 1.0) is similar to loop B, but by 
examining the membrane voltage in Figure |4f), one sees 
that the event termination transitions smoothly from de- 
polarization block back to the resting level without ex- 
hibiting any spikes. This is because the return trip along 
the upper portion of the loop avoids the Hopf bifurcation. 
We clarify this transition in the next section. (Note that, 
for nearby parameter sets, it is possible to observe bursts 
of this type with more spikes at the event onset than are 
shown here.) 

Figure [S] shows a magnification of the upper part of 
Figure [TJ3 that reveals more detail. As the SNIC curve 
continues up from the lower left corner, a codimension- 
two bifurcation known as a saddle-node loop (SL) 
is encountered [111 approximately ([if]o, [Na]i) — 
(14.994, 34.795) mM. At this point, the curve splits into 
two branches corresponding to saddle-node (SN) and ho- 
moclinic bifurcations (HC). The saddle- node branch con- 
tinues to the upper right and forms a cusp with another 



saddle-node branch, while the homoclinic branch contin- 
ues up, curves around, and terminates at a (codimension- 
two) Bogdanov-Takens (BT) point at approximately 
(16.917, 36.714) mM. This point is coincident with the 
upper SN branch, and is also an endpoint of the Hopf 
bifurcation curve (HB). 

Superimposed on this diagram is the upper portion of 
loop D from Figure[3l which represents the termination of 
the burst event. With the increased magnification, it can 
be seen that this part of the ion concentration limit cycle 
indeed does not cross the Hopf curve. Instead, it crosses 
the two saddle-node curves. To clarify the nature of the 
burst termination, we show in the left panel of Figure 
|6] (labeled 'D') the one-dimensional bifurcation diagram 
analogous to Figure [TJ\ for [Na]i = 37.2 mM, along with 
the burst termination portion of the system's trajectory. 
Also included is an inset showing the membrane voltage 
versus time for one complete burst event. It can be seen 
that the termination (arrows) occurs when the trajectory 
tracking the upper stable equilibrium branch encounters 
the second saddle-node bifurcation and drops to the lower 
stable equilibrium branch. 

The additional detail in Figure [5] permits the identi- 
fication of a fifth type of bursting pattern exhibited by 
our system (obtained with p = 0.9, G = 10, e = 0.5, 
kbath = 20, and 7 = 0.25). E labels the upper portion of 
an ion concentration loop that is slightly lower than the D 
loop. It can be seen that this portion of loop E crosses, in 
order, the Hopf, saddle-node, and homoclinic bifurcation 
curves. The remaining portion of loop E (not shown) 
is similar to loop B, and indeed the membrane voltage 
traces are similar, but the termination mechanisms are 
different. Notably, this type of burst terminates after 
a few small-amplitude spikes via homoclinic rather than 
SNIC bifurcation; see Figure [BJ3. 

Finally, we note that we have also identified parameter 
sets that give rise to loops that encircle the BT point in 
Figure [5] while lying entirely above the SL point, as well 
as similar loops occurring upward and to the right of 
this, including ones that straddle only the two SN curves. 
However, the separation of time scales is less clear in 
these cases. And although it might exist, we were not 
able to find a loop that crosses the SNIC and Hopf curves 
and then returns above the saddle-node cusp. Such a 
loop would exhibit a burst termination that transitions 
from depolarization block to the rest state smoothly and 
without a sudden drop in membrane potential. 



IV. DISCUSSION 

In this paper, we have presented a simple model of 
a single neuron with dynamic intra- and extra-cellular 
sodium and potassium concentrations that exhibits pe- 
riodic bursting behavior. Our goal has been to present 
a catalog of the wide variety of qualitatively different 
bursting patterns exhibited by this model based on an 
understanding of the underlying bifurcation structure. 
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FIG. 4: Four qualitatively different bursting patterns corresponding to the four loops shown in Figure [S] In the lower panels, 
solid curves represent [K]o (left vertical axes), and dotted curves represent [Na]i (right vertical axes). The event in (D) recurs 
with a period of approximately 16.5 seconds, but for clarity, only a portion of one such event is shown. 



To do this, we have taken advantage of the large separa- 
tion of time scales inherent in the model: that of spikes, 
which occur on the order of milliseconds, and that of the 
bursting events, which occur on the order of tens of sec- 
onds. We have focused attention on the nature of the 
fast spiking dynamics by freezing the slow ion concentra- 
tion variables and identifying the asymptotic behavior 
of the neuron under these conditions. This clarifies how 
oscillations in the slow variables modulate the neuron's 
excitability and give rise to bursting. 

A crucial aspect of our model that is typically not 
present in comparably simple models (to our knowledge) 
is the inclusion of sodium concentration dynamics. Thus, 
we have two slow variables {[Na]i and [K]o) instead of 
one, and it is therefore possible for the slow system to 
undergo bifurcations to oscillatory states. Analysis com- 
plementary to the work presented here, in which bifurca- 



tions in the slow dynamics are analyzed by removin g th e 
fast dynamics via averaging, has been presented in [18| . 

For this study, we deliberately chose a single neuron 
with the classic Hodgkin-Huxley ionic currents in order 
to focus attention on the role of the ion concentration 
dynamics. Consequently, some quantitative aspects of 
the dynamics (e.g., the value of [K]o for the onset of 
depolarization block in Figure [T]4.) are not directly ap- 
plicable to mammalian systems. More realistic models 
of mammalian systems would require additional currents 
and a network structure with significant synaptic activ- 
ity. (For example, despite the intriguing similarity of 
Figure 4D to phenomena observed in spreading depres- 
sion, it is known that spreading depression requires per- 
sistent sodium channels [13 •) In addition, the absence of 
synaptic background activity means that our model neu- 
ron transitions abruptly from rest directly to bursting as 
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FIG. 5: Magnificaton of the upper portion of Figure[T)3. D and 
E are portions of ion concentration limit cycles as described 
in the text. Abbreviations: SN=saddle-node, BT=Bogdanov- 
Takens, HC=homoclinic, HB=Hopf, SL=Saddle-node loop, 
SNIC=saddle-node infinite-period. 



[K]o increases. In a more realistic setting, we would ex- 
pect to see an intervening regime of tonic spiking due to 
network effects before the onset of bursting. 

Although our approach has concentrated on identifying 
all possible bursting scenarios without regard to physi- 
ological relevance, we note that bursting morphologies 
similar to those that we catalog here have been observed 
in experimental models of epilepsy, including a high- 
potassium model (2I], a high-potassium, low-calcium 
model , and a low- magnesium 4-aminopyridine model 
This suggests that despite the obvious limitations 
of our model, it is capturing aspects of the essential dy- 
namics in these experiments. And, although there has 
been debate about whether a single neuron can undergo 
a "seizure" [H,!!^], our single- neuron results suggest that 
ion concentration dynamics may well play an important 
role in epilepsy. 

It is interesting to note that in order to obtain our 
type D burst, in which the burst termination transitions 
smoothly from depolarization block back to rest without 
exhibiting spikes (Figure [If)), we found it necessary to 
increase the parameter 7, which is inversely proportional 
to the volume of our assumed spherical cell (see the Ap- 
pendix). With a smaller volume, the internal sodium 
concentration is more easily increased to sufficient levels 
to prevent spiking. Thus, our model predicts that this 
type of burst should be most easily observed in smaller 
neurons. Accordingly, this behavior has been seen in hip- 
pocampal interneurons (e.g.. Figure 4a of [1^). Similar 
behavior in other cases, however, may well be due to 
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FIG. 6: (Color) Diagrams clarifying the burst termination 
behavior corresponding to ion concentration loops D and E 
of Figure (5] Insets show complete burst time traces of the 
membrane voltage. In D, [Na]i — 37.2 mM; in E, [Na]i = 
36.6 niM. 



different mechanisms (e.g., calcium [31| or chloride accu- 
mulation [32]), and future work will explore the conse- 
quences of including these and other ions [ssj , as well as 
a more realistic collection of channels. 



V. CONCLUSION 

We have found that a simple model of a single neuron, 
augmented by dynamic intra- and extra-cellular ion con- 
centrations, can display various kinds of periodic burst- 
ing behavior. Our work emphasizes the importance of ion 
concentration homeostasis in the maintenance of the nor- 
mal physiological state, and suggests that a breakdown in 
such homeostatic mechanisms may underlie pathological 
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conditions such as epilepsy. 
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Appendix A: Parameters and Equations for Gating 
Variables 



Parameters not specified above were set as follows: 

Cm — 1 l^Fjcra^ 
QNa = 100 mS/crn^ 
9NaL = 0.0175 mS/crn^ 
qk — 40 mS/crn^ 
9kl = 0.05 mS/cm? 
gciL = 0.05 mSjcrr?. 

The sodium activation was instantaneous, with 
moo {V) -^^^ 



a™(V)+/3„(V) 
^ (Y) - 0-1(^+30) 
'^rnyv ) — l-cxp(-0.1(y+30)) 

The remaining gating variables were governed by 
(j) = 3.0 

a^(V^)=0.07exp(-^) 

Ph {V) — i+cxp(-0.1(y+14)) 
_ 0.01(V+34) 
7 — l-oxp(-0.1(V+34)) 

/3„(V^) = 0.125 exp(-^). 



with 



Appendix B: Derivation of Conversion Factors 

To derive the conversion factor 7, we consider a mem- 
brane current density /, measured in ^A/cm?, and ask: 
how many ions exit a cell through a patch of membrane of 
area A [cm?) in time Ai (sec) due to this current? Recall- 
ing the convention that membrane currents are positive- 
outward, the amount of charge that leaves is 

Ag = {I A) At. 

Assuming monovalent ions, the number of ions is 

-;(^)- 

where q — 1.60 x 10 Coulombs. Hence, the rate at 
which ions pass through the area A is 



dN 
IF 



lA 

q 



Now assume a spherical cell of volume Vin (measured in 
CTO'^) with an intracellular concentration Ci„ of positive 
ions measured in ions/cm^. We have 

- 

where N is the number of ions within the volume. Hence, 



dci, 



dt 



1 

V~ 



dN 
'dt 



1 



lA 



The minus sign has been introduced because a positive 
(outward) current corresponds to a decrease in the num- 
ber of (positive) ions within the cell. 

We wish to express this rate-of-change of concentration 
in mM/msec. The expression on the right-hand side has 
units 



1 



cm^ 



cm'^ 



Coul 



10" 



; ions 1 millimole 
> 

cm^sec msec L 



where Na = 6.02 x 10^'^ is Avogadro's number. Since 1 
millimole / L — 1 mM , we have 



dci, 



dt 



1 



lA 



NAqV^r. 



We write 



with 



dci, 



dt 



1 



7I 



(Bl) 



A 



7 



NAqV„ 



A 



where F = NAq is the Faraday constant. Equation (|B1|) 
is used for the intracellular sodium in Equation 

Now assume that the sphere has radius r, and let A be 
its total surface area. (More correct would be to let A be 
just the area of the channel pores.) Then, 
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_3_ 



Using r — 7 ^m, we obtain 7 — 4.45 x 10^'^. Conversely, 
for 7 = 0.25 and 1.0, we obtain r = 1.24 and 0.31 /im, 
respectively. 

Consider now the concentration of positive ions in the 
extracellular space. The volume in the expression for 
7 must be replaced with the extracellular volume, and 
the negative sign in Equation IBll must be removed. Let 
P = Vin/Vout- Then 7 — 7/?, and we have 



dCn 



dt 



1 

10^ 



7/3/. 



(B2) 



Equation (jB2[) is used for the extracellular potassium in 
Equation (g]). 
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